(*  :Title:	Piecewise Convolution Rule Base  *)

(*  :Authors:	Kevin West and Brian Evans  *)

(*
    :Summary:	To take the continuous-time and discrete-time
		convolution of two piecewise-defined functions and return
		another piecewise-defined function.
 *)

(*  :Context:	SignalProcessing`Support`Convolution`  *)

(*  :PackageVersion:  2.7	*)

(*
    :Copyright:	Copyright 1989, 1990 by
		the Digital Signal Processing Group
		at the Georgia Institute of Technology

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of Georgia Tech or Georgia Institute of
	Technology not be used in advertising or publicity pertaining
	to distribution of the software without specific, written prior
	permission.  Georgia Tech makes no representations about the
	suitability of this software for any purpose.  It is provided
	"as is" without express or implied warranty.
 *)

(*  :History:	*)

(*  :Keywords:	*)

(*  :Source:	{Discrete-Time Signal Processing} by Oppenheim & Schafer  *)

(*  :Warning:	*)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:  This package originally performed piecewise continuous
		  convolution and therefore was part of the analog signal
		  processing packages (spp).  We've added discrete convolution,
		  so it is now in the Support packages.  It's old context
		  was "SignalProcessing`Analog`Convolution`"---  if you
		  ask Mathematica to load Convolution from the analog spp,
		  it will still load.
 *)

(*
    :Functions:	AutoCorrelation
		ConvertFromList
		ConvertToList
		CTPiecewiseConvolution
		DTPiecewiseConvolution
		IntervalQ
		PiecewiseConvolution
		PlotList
		SimplifyList
		ValidIntervalQ
 *)



(*  B E G I N     P A C K A G E  *)

BeginPackage [	"SignalProcessing`Support`Convolution`",
		"SignalProcessing`Support`SigProc`",
		"SignalProcessing`Support`SupCode`",
		If [ TrueQ[$VersionNumber >= 2.0],
		     "Algebra`SymbolicSum`",
		     "Algebra`GosperSum`" ] ]


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ]

			
(*  U S A G E     I N F O R M A T I O N  *)

$ConvolutionDomain::usage = 
	"$ConvolutionDomain is the default domain for all of the routines \
	in the Convolution package (except for DTPiecewiseConvolution and \
	CTPiecewiseConvolution which imply a domain). \
	Use SetConvolutionDomain to reset its value."

Area::usage =
	"Area[ar] is an object representing the area under a Dirac delta \
	function. \
	A Dirac delta function is usually written as C Delta[t - t0], \
	where C is the area under the delta function and t0 is the location \
	of the delta function. \
	In F-interval form, C Delta[t - t0] becomes { Area[C], t0, t0 }."

AutoCorrelation::usage =
	"AutoCorrelation[e, v] computes the autocorrelation of expression e \
	(in the form of a list of F-intervals, an F-interval, or a function) \
	with respect to v. \
	The variable v is treated as either discrete and continuous \
	according to the value of the Domain option. \
	The default domain is the value of $ConvolutionDomain."

ConvertFromList::usage =
	"ConvertFromList[l, t] converts a list of F-intervals in the form \
	{f, t1, t2} to a signal processing expression. \
	The new expression will contain either Delta, CStep, and CPulse \
	functions (continuous-time) or Impulse, Step, and Pulse functions \
	(discrete-time) according to the value of the Domain option. \
	The default domain is the value of $ConvolutionDomain."

ConvertToList::usage =
	"ConvertToList[e, t] converts the expression e to a \
	list of F-intervals in the form {fun, left, right}. \
	The domain of the new expression is established according to the \
	value of the Domain option. \
	The default domain is the value of $ConvolutionDomain."

CTPiecewiseConvolution::usage =
	"CTPiecewiseConvolution[f, g, t] carries out the one-dimensional,
	continuous-time piecewise convolution of f and g with respect to t. \
	See PiecewiseConvolution for valid representation of f and g. \
	Using F-interval notation {function, left_endpoint, right_endpoint}, \
	a Dirac delta is meant when function = Area[<area>] \
	(where <area> equals the area of under the impulse function) and \
	left_endpoint = right_endpoint (which is the location of the impulse)."

DTPiecewiseConvolution::usage = 
	"DTPiecewiseConvolution[f, g, n] carries out the one-dimensional,
	discrete-time piecewise convolution of f and g with respect to n. \
	See PiecewiseConvolution for the valid representations of f and g. \
	Using F-interval notation {function, left_endpoint, right_endpoint}, \
	an impulse (Kronecker delta) function is meant when function = <value> \
	(where <value> equals the strength of the impulse function) and \
	left_endpoint = right_endpoint (which is the location of the impulse)."

IntervalQ::usage =
	"IntervalQ[x] returns True if the argument x is an interval, and \
	gives False otherwise. \
	An interval has the form {f, t1, t2} such that none of f, t1, or t2 \
	are lists. \
	See also ValidIntervalQ."

PiecewiseConvolution::usage =
	"PiecewiseConvolution[f, g, v] convolves the piecewise functions \
	f and g with respect to the variable v. \
	In this context, a function is represented in a piecewise fashion: \
	(1) as F-intervals of the form {fun, left, right}, \
	(2) as a list of F-intervals, or (3) as an expression. \
	An F-interval has the form {function, left_endpoint, right_endpoint}. \
	The F-interval notation represents a finite-extent function or \
	sequence when the endpoints do not equal infinity. \
	The result of convolution is always returned as a list of F-intervals. \
	The variable v is treated as either discrete and continuous \
	according to the value of the Domain option. \
	The default domain is the value of $ConvolutionDomain, \
	which can be reset by SetConvolutionDomain. \
	See also CTPiecewiseConvolution and DTPiecewiseConvolution."

PlotList::usage =
	"PlotList[l, {v, vstart, vend}] plots the piecewise function l vs. v. \
	The variable v is treated as either discrete and continuous \
	according to the value of the Domain option. \
	The default domain is the value of $ConvolutionDomain."

SetConvolutionDomain::usage =
	"SetConvolutionDomain[domain] establishes the default convolution \
	domain as either Continuous or Discrete. \
	It returns domain on success, and Null on error."

SimplifyList::usage =
	"SimplifyList[interval_list, v] takes a list of F-intervals \
	each having the form {f, t1, t2} and returns a simplified sorted \
	list of F-intervals. \
	The interval_list is treated as either a continuous-time or a \
	discrete-time function according the value of the Domain option. \
	The variable v is treated as either discrete and continuous \
	according to the value of the Domain option. \
	The default domain is the value of $ConvolutionDomain."

ValidIntervalQ::usage =
	"ValidIntervalQ[i] returns True if i is a valid F-interval, and \
	gives False otherwise. \
	An F-interval is an interval of the form {f, leftendp, rightendp} \
	where f is a function and leftendp < rightendp unless f is an \
	impulse function in which case leftendp must equal rightendp. \
	In a valid F-interval, leftendp and rightenp must be numbers. \
	See also IntervalQ."

(*  E N D     U S A G E     I N F O R M A T I O N  *)



Begin[ "`Private`" ]


(*  S Y S T E M     S E T T I N G S  *)

$ConvolutionDomain = Continuous

Protect[ $ConvolutionDomain ]


(*  M E S S A G E S  *)

AutoCorrelation::domain =
	"The option Domain must be equal to Discrete or Continuous."
ConvertFromList::domain =
	"The option Domain must be equal to Discrete or Continuous."
ConvertToList::domain =
	"The option Domain must be equal to Discrete or Continuous."
CTPiecewiseConvolution::argct =
	"CTPiecewiseConvolution requires exactly three arguments."
CTPiecewiseConvolution::badint =
	"The interval `` is not valid and has been deleted."
DTPiecewiseConvolution::argct =
	"DTPiecewiseConvolution requires exactly three arguments."
DTPlotList::badargs =
	"DTPlotList requires the format DTPlotList[list,{n,n1,n2}]"
PiecewiseConvolution::argct =
	"PiecewiseConvolution requires at least three arguments."
PiecewiseConvolution::domain =
	"The option Domain must be equal to Discrete or Continuous."
PlotList::badargs =
	"PlotList requires the format PlotList[list,{t,t1,t2}]"
PlotList::domain =
	"The option Domain must be equal to Discrete or Continuous."
SimplifyList::domain =
	"The option Domain must be equal to Discrete or Continuous."
SetConvolutionDomain::domain =
	"The option Domain must be equal to Discrete or Continuous."
ValidIntervalQ::domain =
	"The option Domain must be equal to Discrete or Continuous."

(*  E N D     M E S S A G E S  *)


(*  S U P P O R T I N G     R O U T I N E S  *)


(*  AreaOfImpulse -- return the impulse area at a particular point *)
AreaOfImpulse[list_, a_] :=
   Block[ {areaofimp = 0, i},
          Do[ If [ Head[list[[i]][[1]]] == Area,
                   areaofimp += list[[i]][[1]][[1]] * ImpulseHere[list[[i]], a],
                   Null ],
              {i, Length[list]} ];
          areaofimp ]

(*  CallFunction  *)
SetAttributes[CallFunction, {HoldFirst}]

CallFunction[ message_, discretefun_, continuousfun_, oplist_List, args___ ] :=
	Block [	{domain},
		domain = Replace[Domain, oplist];
		Which [	SameQ[domain, Discrete],
			  Apply[discretefun, {args}],
			SameQ[domain, Continuous],
			  Apply[continuousfun, {args}],
			True,
			  message ] ]


(*  ConvertToList routines  *)

Options[ConvertToList] := { Domain :> $ConvolutionDomain }

ConvertToList[f_, v_, options___] :=
	CallFunction[ Message[ConvertToList::domain],
		      DTConvertToList, CTConvertToList,
		      ToList[options] ~Join~ Options[ConvertToList],
		      f, v ]

CTConvertToList[f_, t_] := CTSimplifyList[ UnsimpCTL[f, t], t ]
DTConvertToList[f_, n_] := DTSimplifyList[ UnsimpCTL[f, n], n ]

UnsimpCTL[f_?IntervalQ, t_] := {f}
UnsimpCTL[f_List, t_] := f
UnsimpCTL[f_ + g_, t_] := Join[UnsimpCTL[f,t], UnsimpCTL[g,t]]
UnsimpCTL[0, t_] := {}
UnsimpCTL[Delta[t_ + a_.] ar_., t_] :=
	{{Area[ Limit[ar, t -> -a] ], -a, -a}} /; FreeQ[a,t]
UnsimpCTL[Delta[-t_ + a_.] ar_., t_] :=
	{{Area[ Limit[ar, t -> a] ], a, a}} /; FreeQ[a,t]
UnsimpCTL[Impulse[t_ + a_.] ar_., t_] :=
	{{Limit[ ar, t -> -a], -a, -a}} /; FreeQ[a,t]
UnsimpCTL[Impulse[-t_ + a_.] ar_., t_] :=
	{{Limit[ar, t -> a], a, a}} /; FreeQ[a,t]
UnsimpCTL[CStep[t_ + a_.] f_., t_] := {{f,-a,Infinity}} /; FreeQ[a, t]
UnsimpCTL[CStep[-t_ + a_.] f_., t_] := {{f,-Infinity,a}} /; FreeQ[a, t]
UnsimpCTL[Step[t_ + a_.] f_., t_] := {{f,-a,Infinity}} /; FreeQ[a, t]
UnsimpCTL[Step[-t_ + a_.] f_., t_] := {{f,-Infinity,a}} /; FreeQ[a, t]
UnsimpCTL[CPulse[l_, t_ + a_.] f_.,t_] := {{f, -a, l - a}} /; FreeQ[{l,a}, t]
UnsimpCTL[CPulse[l_, -t_ + a_.] f_.,t_] := {{f, a - l, a}} /; FreeQ[{l,a}, t]
UnsimpCTL[Pulse[l_, t_ + a_.] f_.,t_] := {{f, -a, l - 1 - a}} /; FreeQ[{l,a}, t]
UnsimpCTL[Pulse[l_, -t_ + a_.] f_.,t_] := {{f, a - l + 1, a}} /; FreeQ[{l,a}, t]
UnsimpCTL[f_ term_Plus, t_] := UnsimpCTL[ Distribute[f term], t ] /;
	! FreeQ[term, t] && ! SameQ[f term, Distribute[f term]]
UnsimpCTL[f_, t_] := {{f,-Infinity,Infinity}}

(*  ConvertFromList routines  *)

Options[ConvertFromList] := { Domain :> $ConvolutionDomain }

ConvertFromList[list_, v_, options___] :=
	CallFunction[ Message[ConvertFromList::domain],
		      DTConvertFromList, CTConvertFromList,
		      ToList[options] ~Join~ Options[ConvertFromList],
		      list, v ]

CTConvertFromList[list_, t_] :=
	Apply[Plus, Map[CTConvertFromInterval[#1, t]&, list]]
DTConvertFromList[list_, t_] :=
        Apply[Plus, Map[DTConvertFromInterval[#1, t]&, list]]

CTConvertFromInterval[{Area[ar_], a_, a_}, t_] := ar Delta[t - a]
CTConvertFromInterval[{f_, -Infinity, Infinity}, t_] := f
CTConvertFromInterval[{f_, a_, Infinity}, t_] := f CStep[t - a]
CTConvertFromInterval[{f_, -Infinity, a_}, t_] := f CStep[a - t]
CTConvertFromInterval[{f_, a_, b_}, t_] := f CPulse[b - a, t - a]

DTConvertFromInterval[{Area[ar_], a_, a_}, t_] := ar Impulse[t-a]  
DTConvertFromInterval[{ar_, a_, a_}, t_] := ar Impulse[t-a]
DTConvertFromInterval[{f_, -Infinity, Infinity}, t_] := f
DTConvertFromInterval[{f_, a_, Infinity}, t_] := f Step[t - a]
DTConvertFromInterval[{f_, -Infinity, a_}, t_] := f Step[a - t]
DTConvertFromInterval[{f_, a_, b_}, t_] := f Pulse[b - a + 1, t - a]


(*  ConvolveTwoIntervals  *)
ConvolveTwoIntervals[pair_List, t_, time_] :=
	ConvolveTwoIntervals[ pair[[1]], pair[[2]], t, time]

ConvolveTwoIntervals[{Area[ar1_], a_, a_}, {Area[ar2_], b_, b_}, t_, time_] :=
	{{Area[ar1 ar2], a+b, a+b}} 
ConvolveTwoIntervals[{Area[ar_], a_, a_}, {f_, b_, c_}, t_, time_] :=
	{{ar (f /. t -> t-a), b+a, c+a}}
ConvolveTwoIntervals[{f_, a_, b_}, {Area[ar_], c_, c_}, t_, time_] :=
	{{ar (f /. t -> t-c), a+c, b+c}}
ConvolveTwoIntervals[{f_, -Infinity, t2_}, {g_, t3_, Infinity}, t_, time_] :=
	ConvolveTwoIntervals[{g, t3, Infinity}, {f, -Infinity, t2}, t, time] /;
	t3 > -Infinity
	  (* Jan-92 JMc: make sure left-sided signal is g(t) *)
ConvolveTwoIntervals[{f_, t1_, t2_}, {g_, t3_, t4_}, t_, time_] :=
	ConvolveTwoIntervals[{g, t3, t4}, {f, t1, t2}, t, time] /;
	(t4 - t3) > (t2 - t1)

ConvolveTwoIntervals[{f_, t1_, t2_}, {g_, t3_, t4_}, t_, Continuous] := 
	Block [	{convlist, f2, g2, T, ta, tb, tc, td},
		f2 = f /. t -> T;
		g2 = g /. t -> t - T;
		ta = If [ t1 > -Infinity, t1+t4,  t1 ];
		tb = If [ t2 <  Infinity, t2+t3,  t2 ];
		tc = If [ t3 > -Infinity, t -t3, -t3 ];
		td = If [ t4 <  Infinity, t -t4, -t4 ];
		convlist = {};
		If [ ta > -Infinity,
		     AppendTo[convlist,
		              {Integrate[f2 g2, {T,t1,tc}], t1+t3, ta}] ];
		If [ ta < Infinity && tb > -Infinity,
		     AppendTo[convlist,
		              {Integrate[f2 g2, {T,td,tc}], ta, tb}] ];
		If [ tb < Infinity,
		     AppendTo[convlist,
		              {Integrate[f2 g2, {T,td,t2}], tb, t2+t4}] ];
		convlist ]

ConvolveTwoIntervals[{f_, t1_, t2_}, {g_, t3_, t4_}, t_, Discrete] := 
	Block [	{convlist, f2, g2, sumhead, summation, tt1, tt2, T, ta, tb, tc, td},
		f2 = f /. t -> T;
		g2 = g /. t -> t - T;
		ta = If [ t1 > -Infinity, t1+t4,  t1 ];
		tb = If [ t2 <  Infinity, t2+t3,  t2 ];
		tc = If [ t3 > -Infinity, t -t3, -t3 ];
		td = If [ t4 <  Infinity, t -t4, -t4 ];
		sumhead = If [ TrueQ[$VersionNumber >= 2.0],
			       SymbolicSum,
			       GosperSum ];
		summation = Cancel[sumhead[f2 g2, {T, tt1, tt2}]];
		convlist = {};
		If [ ta > -Infinity,
		     AppendTo[convlist,
		              {summation /. tt1 -> t1 /. tt2 -> tc, t1+t3, ta}] ];
		If [ ta < Infinity && tb > -Infinity,
		     AppendTo[convlist,
		              {summation /. tt1 -> td /. tt2 -> tc, ta+1, tb}] ];
		If [ tb < Infinity,
		     AppendTo[convlist,
		              {summation /. tt1 -> td /. tt2 -> t2, tb+1, t2+t4+1}] ];
		convlist ]

(*  ConvTwoLists  *)
PairIntervalAndList[int_, list_] := Map[ {int, #}&, list ]

CTConvTwoLists[a_, b_, v_] := ConvTwoLists[a, b, v, Continuous]
DTConvTwoLists[a_, b_, v_] := ConvTwoLists[a, b, v, Discrete]

ConvTwoLists[{}, b_, v_, time_] := {}
ConvTwoLists[a_, {}, v_, time_] := {}
ConvTwoLists[a_, b_, v_, time_] :=
	Block [	{convolution, groupedinto3s, nestedlist, pairedlists},
		nestedlist = Flatten[ Map[ PairIntervalAndList[#, b]&, a ] ];
		pairedlists = Partition[ Partition[nestedlist, 3], 2 ];
		convolution = Map[ ConvolveTwoIntervals[#, v, time]&,
				   pairedlists ];
		groupedinto3s = Partition[ Flatten[convolution], 3 ];
		SimplifyList[ RemoveBadInterval[groupedinto3s, time], v,
			      Domain -> time ] ]

(*  ImpulseHere -- return 1 if the points coincide *)
ImpulseHere[{Area[ar_], a_, a_}, a_] := 1
ImpulseHere[___] := 0

(*  InIntervalQ; InInterval returns a 1 if the intervals cross  *)    
InIntervalQ[{f_, left_, right_}, a_, b_] :=
	TrueQ[ N[MyLessEqual[left, a, b, right]] ]
InInterval[args__] := If [ InIntervalQ[args], 1, 0 ]

(*  IntervalQ  *)
IntervalQ[{f_, a_, b_}] := ! ( ListQ[f] || ListQ[a] || ListQ[b] )
IntervalQ[___] := False

(*  ListOfEndpoints  *)
getsecthird[ {a_, b_, c_} ] := ToCollection[b, c]

ListOfEndpoints[list_List] := Sort[ Union[ Map[getsecthird, list] ], MyLess ]


(*  MyLess  *)
MyLess[a_, a_] := False	
MyLess[-Infinity, Infinity] := True
MyLess[Infinity, -Infinity] := False
MyLess[Infinity, a_] := False
MyLess[a_, Infinity] := True
MyLess[a_, -Infinity] := False
MyLess[-Infinity, a_] := True
MyLess[a_, b_] := N[a] < N[b]
MyLess[a_, b_, rest__] := MyLess[a, b] && MyLess[b, rest]

(*  MyLessEqual  *)
MyLessEqual[a_, a_] := True
MyLessEqual[a_, b_] := MyLess[a, b]
MyLessEqual[a_, b_, rest__] := MyLessEqual[a, b] && MyLessEqual[b, rest]

(*  RemoveBadInterval  *)
RemoveBadInterval[oldlist_, Discrete] := Select[oldlist, DTValidIntervalQ]
RemoveBadInterval[oldlist_, Continuous] := Select[oldlist, CTValidIntervalQ]

(*  RemoveDomain  *)
keepq[x_] := SameQ[Domain, Replace[Domain, x]]
RemoveDomain[oplist_] := Select[oplist, keepq]

(*  ValueOfInterval -- return the value of a particular interval  *)
ValueOfInterval[list_, a_, b_] :=
	Apply[Plus, Map[First, Select[list, InIntervalQ[#1, a, b]&]]]

(*  ValidIntervalQ  *)
Options[ValidIntervalQ] := { Domain :> $ConvolutionDomain }

ValidIntervalQ[list_, options___] := 
	CallFunction[ Message[ValidIntervalQ::domain],
		      DTValidIntervalQ, CTValidIntervalQ,
		      ToList[options] ~Join~ Options[ValidIntervalQ],
		      list ]

CTValidIntervalQ[i_] := False /; Not[IntervalQ[i]]
CTValidIntervalQ[{Area[ar_], a_, b_}] := TrueQ[ (a == b) && (ar != 0) ]
CTValidIntervalQ[{f_, a_, b_}] := TrueQ[ N[(b > a) && ! SameQ[f, 0]] ]

DTValidIntervalQ[i_] := False /; Not[IntervalQ[i]]
DTValidIntervalQ[{Area[ar_], a_, b_}] :=
	SameQ[a, b] && (! InfinityQ[b]) && (! SameQ[ar, 0])
DTValidIntervalQ[{f_, a_, a_}] := (! InfinityQ[a]) && (! SameQ[f, 0])
DTValidIntervalQ[{f_, a_, b_}] := TrueQ[ N[(b >= a) && (! SameQ[f, 0])] ]


(*  C O N V O L U T I O N     R O U T I N E S  *)


(*  AutoCorrelation  *)

Options[AutoCorrelation] := { Domain :> $ConvolutionDomain }

AutoCorrelation[f_, v_, options___] :=
	CallFunction[ Message[AutoCorrelation::domain],
		      DTPiecewiseConvolution, CTPiecewiseConvolution,
		      ToList[options] ~Join~ Options[AutoCorrelation],
		      f, FlipInTime[f, v], v ]

FlipInTime[int_, v_] := FlipInterval[int, v] /; IntervalQ[int]
FlipInTime[list_List, v_] := Reverse[ Map[FlipInterval[#1, v]&, list] ]
FlipInTime[f_, v_] := f /. v -> -v

FlipInterval[{f_, v1_, v2_}, v_] := {(f /. v -> -v), -v2, -v1}

(*  PiecewiseConvolution  *)
Options[PiecewiseConvolution] := { Domain :> $ConvolutionDomain }

PiecewiseConvolution[f_, g_, v_, options___] :=
	CallFunction[ Message[PiecewiseConvolution::domain],
		      DTPiecewiseConvolution, CTPiecewiseConvolution,
		      ToList[options] ~Join~ Options[PiecewiseConvolution],
		      f, g, v ]

PiecewiseConvolution[___] := Message[PiecewiseConvolution::argct]

Unprotect[CConvolve, Convolve]

(*  CTPiecewiseConvolution  *)
CConvolve/: TheFunction[ CConvolve[t_Symbol][x_, y_] ] :=
	CTPiecewiseConvolution[x, y, t]
CConvolve/: TheFunction[ CConvolve[t_Symbol][x_, y_, rest__] ] :=
	TheFunction[ CConvolve[t][ CTPiecewiseConvolution[x, y, t], rest ] ]

CTPiecewiseConvolution[f_, g_, t_] :=
	CTConvTwoLists[ CTConvertToList[f, t], CTConvertToList[g, t], t ]
CTPiecewiseConvolution[___] := Message[CTPiecewiseConvolution::argct]

(*  DTPiecewiseConvolution  *)
Convolve/: TheFunction[ Convolve[n_Symbol][x_, y_] ] :=
	DTPiecewiseConvolution[x, y, n]
Convolve/: TheFunction[ Convolve[n_Symbol][x_, y_, rest__] ] :=
	TheFunction[ Convolve[n][ DTPiecewiseConvolution[x, y, n], rest ] ]

DTPiecewiseConvolution[f_, g_, n_] :=
	DTConvTwoLists[ DTConvertToList[f, n], DTConvertToList[g, n], n ]
DTPiecewiseConvolution[___] := Message[DTPiecewiseConvolution::argct]

Protect[CConvolve, Convolve]

(*  P L O T T I N G     R O U T I N E S  *)

(*  PlotList  *)
Options[PlotList] := { Domain :> $ConvolutionDomain }

PlotList[list_, {v_, v1_, v2_}, options___] :=
	Block [	{oplist},
		oplist = ToList[options] ~Join~ Options[PlotList];
		CallFunction[ Message[PlotList::domain],
			      DTPlotList, CTPlotList, oplist,
			      list, {v, v1, v2}, RemoveDomain[oplist] ] ]

(* CTPlotList  *)
CTPlotList[int_, {t_, t1_, t2_}, options___] :=
        PlotList[{int}, {t, t1, t2}, options] /; IntervalQ[int]
CTPlotList[list_List, {t_, t1_, t2_}, options___] :=
	SignalPlot[CTConvertFromList[list, t], {t, t1, t2}, options]
CTPlotList[fun_, {t_, t1_, t2_}, options___ ] :=
	SignalPlot[fun, {t, t1, t2}, options]

(* DTPlotList *)
DTPlotList[int_, {n_, n1_, n2_}, options___] :=
	DTPlotList[{int}, {n, n1, n2}, options] /; IntervalQ[int]
DTPlotList[list_List, {n_, n1_, n2_}, options___] :=
	SequencePlot[DTConvertFromList[list, n], {n, n1, n2}, options]
DTPlotList[fun_, {n_, n1_, n2_}, options___] :=
	SequencePlot[fun, {n, n1, n2}, options]


(*  S I M P L I F I C A T I O N     R O U T I N E S  *)

(*  SimplifyList  *)
Options[SimplifyList] := { Domain :> $ConvolutionDomain }

SimplifyList[list_, v_, options___] :=
	CallFunction[ Message[SimplifyList::domain],
		      DTSimplifyList, CTSimplifyList,
		      ToList[options] ~Join~ Options[SimplifyList],
		      list, v ]

(*  CTSimplifyList  *)
CTSimplifyList[{}, t_] := {}
CTSimplifyList[list_, t_] :=
	Block[ { curpoint, k, nextpoint, points, simplelist},
		points = ListOfEndpoints[list];		
		nextpoint = points[[1]];
		simplelist = {{ Area[AreaOfImpulse[list, nextpoint]],
				nextpoint, nextpoint }};
		Do [ curpoint = nextpoint;
		     nextpoint = points[[k]];
		     AppendTo[ simplelist,
			       { ValueOfInterval[list, curpoint, nextpoint],
				 curpoint, nextpoint } ];
		     AppendTo[ simplelist,
			       { Area[AreaOfImpulse[list, nextpoint]],
	       	                 nextpoint, nextpoint } ],
		     {k, 2, Length[points]} ];
		RemoveBadInterval[SPSimplify[simplelist, Variables -> {t}],
				  Continuous] ]

(*  DTSimplifyList  *)
OpenList[{Area[A_], a_, a_}] := {A, a, a+1}
OpenList[{Area[A_], a_, b_}] := {} /; ! SameQ[a, b]
OpenList[{f_, a_, b_}] := {f, a, b+1}
CloseList[{f_, a_, b_}] := {f, a, b-1}

DTSimplifyList[{}, n_] := {}
DTSimplifyList[list_, n_] :=
	Block [	{curpoint, ival, k, lastpoint, nextpoint, numpoints,
		 openlist, points, simplelist = {}, simplified},

		openlist = Map[OpenList, list];
		points = ListOfEndpoints[openlist];
		numpoints = Length[points];
		nextpoint = points[[1]];
		Do [ curpoint = nextpoint;
		     nextpoint = points[[k]];
		     ival = ValueOfInterval[openlist, curpoint, nextpoint];
		     AppendTo[ simplelist, {ival, curpoint, nextpoint} ],
		     {k, 2, numpoints} ];
		simplified = SPSimplify[Map[CloseList, simplelist],
					Variables -> {n}];
		simplelist = RemoveBadInterval[simplified, Discrete];

		openlist = {};
		While [	Length[simplelist] > 0,
			curpoint = First[ simplelist ];
			simplelist = Rest[ simplelist ];

			If [ ! SameQ[ curpoint[[2]], curpoint[[3]] ],
			     AppendTo[ openlist, curpoint ];
			     curpoint = {},

			     If [ Length[openlist] > 0,
				  lastpoint = Last[openlist];
				  If [ SameQ[curpoint[[2]], 1+lastpoint[[3]]] &&
				       SameQ[curpoint[[1]] /. n->curpoint[[2]],
				           lastpoint[[1]] /. n->curpoint[[2]]],
				       openlist = Append[ Drop[openlist, -1],
							  { lastpoint[[1]],
							    lastpoint[[2]],
							    curpoint[[2]] } ];
				       curpoint = {} ] ];

			     If [ Length[simplelist] > 0 &&
				    Length[curpoint] > 0,
				  nextpoint = First[simplelist];
				  If [ SameQ[curpoint[[2]], nextpoint[[2]]-1] &&
				       SameQ[curpoint[[1]] /. n->curpoint[[2]],
				          nextpoint[[1]] /. n->curpoint[[2]]],
				       AppendTo[ openlist,
						 { nextpoint[[1]],
						   nextpoint[[2]] - 1,
						   nextpoint[[3]] } ];
				       simplelist = Rest[simplelist];
				       curpoint = {} ] ];

			     If [ Length[curpoint] > 0,
				  AppendTo[openlist, curpoint] ] ] ];

		lastpoint = Last[openlist];
		While [ SameQ[First[lastpoint] /. n -> lastpoint[[3]], 0] &&
			! SameQ[lastpoint[[3]], Infinity],
			openlist = If [ SameQ[lastpoint[[2]], lastpoint[[3]]],
					Drop[openlist, -1],
		    			Append[ Drop[openlist, -1],
						{ lastpoint[[1]],
						  lastpoint[[2]],
						  lastpoint[[3]] - 1 } ] ];
			lastpoint = Last[ openlist ] ];

		openlist ]

(*  SetConvolutionDomain  *)
SetConvolutionDomain[domain_] :=
	Block [	{},
		Unprotect[ $ConvolutionDomain ];
		$ConvolutionDomain = domain;
		Protect[ $ConvolutionDomain ];
		domain ] /;
	SameQ[domain, Discrete] || SameQ[domain, Continuous]

SetConvolutionDomain[bogusdomain_] := Message[SetConvolutionDomain::domain]


(*  E N D     P A C K A G E  *) 

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell1 ];
     On[ General::spell ] ]


(*  A L I A S E S  *)

PiecewisePlot = PlotList
PiecewisePlot::usage = PlotList::usage


(*  H E L P     I N F O R M A T I O N  *)

Block [	{newfuns},
	 newfuns =
	 { AutoCorrelation,		ConvertFromList,	ConvertToList,
	   CTPiecewiseConvolution,	DTPiecewiseConvolution, IntervalQ,
	   PiecewiseConvolution,	PiecewisePlot,		PlotList,
	   SimplifyList,		ValidIntervalQ };
	 Combine[ SPfunctions, newfuns ];
	 Apply[ Protect, newfuns ] ]


(*  E N D I N G     M E S S A G E  *)

Print["Piecewise convolution rules have been loaded."]
Null
